Home Range Analysis

Home Range Analysis

The data presented here was sampled from movebank.org and was originally published here:

Carrasco-Harris MF, Bowman D, Reichling S, Cole JA. 2020. Spatial ecology of copperhead snakes (Agkistrodon contortrix) in response to urban park trails. J Urban Ecol. 6(1):juaa007. https://doi.org/10.1093/jue/juaa007

In their study, the authors measured the spatial activity of copperheads (Agkistrodon contortrix) around pedestrian trails in an urban forest (Overton Park, Memphis, TN) using radio telemetry. They found that sex and season, not distance to the nearest trail, affected the distance that snakes moved and concluded that copperheads may be tolerant of low-level human disturbances.

\(~\)

1.1 Bring in data

\(~\)

Import Dataset with read.csv

snake.data <- read.csv("data/snakes.csv")
snake <- snake.data %>%
    filter(individual.local.identifier == c("ID_4241","ID_4270","ID_7270","ID_4266"))

## Warning in individual.local.identifier == c("ID_4241", "ID_4270", "ID_7270", :
## longer object length is not a multiple of shorter object length

\(~\)

I decided to subset for just 4 individuals that I thought gave a good representation of the data. These snakes each primarily resided in different sections of the park and displayed variation in the sizes of their respective home ranges.

\(~\)

head(snake)

##      event.id visible               timestamp location.long location.lat
## 1 13872701454    true 2017-05-16 15:26:00.000     -89.98502     35.14397
## 2 13872701458    true 2017-05-29 16:00:00.000     -89.98503     35.14245
## 3 13872701462    true 2017-06-15 14:49:00.000     -89.98473     35.14280
## 4 13872701466    true 2017-07-06 14:57:00.000     -89.98478     35.14288
## 5 13872701470    true 2017-07-20 14:15:00.000     -89.98483     35.14258
## 6 13872701474    true 2017-08-07 14:35:00.000     -89.98482     35.14283
##   study.specific.measurement       sensor.type individual.taxon.canonical.name
## 1      2.67546989:29.6247595 radio-transmitter          Agkistrodon contortrix
## 2    1.244122543:37.94730962 radio-transmitter          Agkistrodon contortrix
## 3    11.34398646:23.69300256 radio-transmitter          Agkistrodon contortrix
## 4    1.655033961:26.28595886 radio-transmitter          Agkistrodon contortrix
## 5    5.386954834:29.09363462 radio-transmitter          Agkistrodon contortrix
## 6     3.28874227:18.73814023 radio-transmitter          Agkistrodon contortrix
##   tag.local.identifier individual.local.identifier
## 1              ID_4241                     ID_4241
## 2              ID_4241                     ID_4241
## 3              ID_4241                     ID_4241
## 4              ID_4241                     ID_4241
## 5              ID_4241                     ID_4241
## 6              ID_4241                     ID_4241
##                                              study.name utm.easting
## 1 Spatial ecology of copperhead snakes in an urban park    228048.4
## 2 Spatial ecology of copperhead snakes in an urban park    228041.8
## 3 Spatial ecology of copperhead snakes in an urban park    228070.3
## 4 Spatial ecology of copperhead snakes in an urban park    228066.1
## 5 Spatial ecology of copperhead snakes in an urban park    228060.5
## 6 Spatial ecology of copperhead snakes in an urban park    228062.9
##   utm.northing utm.zone
## 1      3893089      16N
## 2      3892920      16N
## 3      3892958      16N
## 4      3892968      16N
## 5      3892935      16N
## 6      3892962      16N

\(~\)

1.2 Looking for Outliers

qaqc_plot <- ggplot() + geom_point(data=snake, 
                                   aes(utm.easting,utm.northing,
                                       color=individual.local.identifier)) +
                        labs(x="Easting", y="Northing") +
                        guides(color=guide_legend("Identifier"))

ggplotly(qaqc_plot)

\(~\)

The points appear to be have a similar spread for each individual without any clear outliers among them (at least to my eyes).

\(~\)

lapply(split(snake, snake$individual.local.identifier), 
       function(x)write.csv(x, file = paste(x$individual.local.identifier[1],".csv", sep = ""), row.names = FALSE))

## $ID_4241
## NULL
## 
## $ID_4266
## NULL
## 
## $ID_4270
## NULL
## 
## $ID_7270
## NULL

files <- list.files(path = ".", pattern = "[ID_]+[3413-7274]+", full.names = TRUE)

\(~\)

1.3 Analysis

utm_points <- cbind(snake$utm.easting, snake$utm.northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=16 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster2 <- openmap(c(max(proj_lat.lon$y)+0.003, min(proj_lat.lon$x)-0.003), 
                  c(min(proj_lat.lon$y)-0.003, max(proj_lat.lon$x)+0.003), 
                  type = "bing")
raster <- openmap(c(max(proj_lat.lon$y)+0.01, min(proj_lat.lon$x)-0.01), 
                  c(min(proj_lat.lon$y)-0.01, max(proj_lat.lon$x)+0.01), 
                  type = "bing")
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")
raster2_utm <- openproj(raster2, 
              projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")

\(~\)

I decided to make a second raster that still takes into account the location of the points when choosing the maximum and minimum coordinates of the resulting map, but with a smaller window to get us slightly more zoom over our park (which represents such a small area).

Here’s the map made with the raster inputs from Dr. Gentry’s code:

\(~\)

autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="right") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=snake, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 1, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

\(~\)

And with a constrained window:

\(~\)

autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="right") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=snake, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 1, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

\(~\)

1.4 Home Range Analysis

\(~\)

1.4.1 Minimum Convex Polygon

\(~\)

For some reason I was having an issue with the permit status (?) using the pbapply package. For whatever reason and completely mysteriously to me, the first three lines of code in this chunk fixed the issue for me

\(~\)

library(maptools)
if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit() #This resolved the following error message from pbapply: "Error: isTRUE(gpclibPermitStatus()) is not TRUE"

## [1] TRUE

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual.local.identifier)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), alpha=0.5, fill = "red") +
    geom_point(data=mcp.points, aes(x=x, y=y)) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

\(~\)

1.4.2 Kernel_Density Estimation

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 95)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual.local.identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.5, fill = "red") +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

\(~\)

1.4.3 Brownian Bridge Movement

snake_4241 <- read.csv("ID_4241.csv")
date <- as.POSIXct(strptime(as.character(snake_4241$timestamp),"%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"))
snake_4241$date <- date
snake_4241.reloc <- cbind.data.frame(snake_4241$utm.easting, snake_4241$utm.northing,
                                as.vector(snake_4241$individual.local.identifier),
                                as.POSIXct(date))
colnames(snake_4241.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(snake_4241.reloc, date=date, id="snake_4241")
sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
snake_4241.traj <- kernelbb(trajectory, sig1 = .7908, sig2 = 58, grid = 100)
bb_ver <- getverticeshr(snake_4241.traj, 58)
bb_poly <- fortify(bb_ver, region = "id", 
                   proj4string = CRS("+proj=utm +zone=16+
                                     datum=WGS84 +units=m +no_defs"))
colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")
bb_image <- crop(snake_4241.traj, bb_ver, 
                 proj4string = CRS("+proj=utm +zone=16 +
                                   datum=WGS84 +units=m +no_defs"))
bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85,  y=0.95,
                   gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
bb.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_tile(data=bb_image, 
            aes(x=bb_image@coords[,1], y=bb_image@coords[,2],
            fill = bb_image@data$ud)) +
  geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "black", fill = NA) +
  scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
  labs(x="Easting (m)", y="Northing (m)", title="4241") +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
bb.plot

\(~\)

snake_4270 <- read.csv("ID_4270.csv")
date <- as.POSIXct(strptime(as.character(snake_4270$timestamp),"%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"))
snake_4270$date <- date
snake_4270.reloc <- cbind.data.frame(snake_4270$utm.easting, snake_4270$utm.northing,
                                as.vector(snake_4270$individual.local.identifier),
                                as.POSIXct(date))
colnames(snake_4270.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(snake_4270.reloc, date=date, id="snake_4270")
sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
snake_4270.traj <- kernelbb(trajectory, sig1 = .7908, sig2 = 58, grid = 100)
bb_ver <- getverticeshr(snake_4270.traj, 24.5)
bb_poly <- fortify(bb_ver, region = "id", 
                   proj4string = CRS("+proj=utm +zone=16+
                                     datum=WGS84 +units=m +no_defs"))
colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")
bb_image <- crop(snake_4270.traj, bb_ver, 
                 proj4string = CRS("+proj=utm +zone=16 +
                                   datum=WGS84 +units=m +no_defs"))
bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85,  y=0.95,
                   gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
bb.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_tile(data=bb_image, 
            aes(x=bb_image@coords[,1], y=bb_image@coords[,2],
            fill = bb_image@data$ud)) +
  geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "black", fill = NA) +
  scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
  labs(x="Easting (m)", y="Northing (m)", title="4270") +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
bb.plot

\(~\)

Interestingly, even though 4270 had a home range that was fairly comparable to 4241 in size (at least by mcp), adding the timestamp data and generating this ‘heat map’ of activity shows that 4270’s activity was generally concentrated in a much denser area.

\(~\)